Imports System.Windows.Forms
Imports System.Threading
Imports MatrixDataStructures

Public Class InhomogenityCorrection
    Sub New()

        ' This call is required by the Windows Form Designer.
        InitializeComponent()

        ' Add any initialization after the InitializeComponent() call.

    End Sub
    Sub New(ByVal Data As MatrixDataStructures.Matrix3DInt16)

        ' This call is required by the Windows Form Designer.
        InitializeComponent()
        Me.ShowIcon = False
        ' Add any initialization after the InitializeComponent() call.
        Me.DataInt = New MatrixDataStructures.Matrix3DInt16(Data.x, Data.y, Data.z)
        For i As Integer = 0 To Data.x - 1
            For j As Integer = 0 To Data.y - 1
                For k As Integer = 0 To Data.z - 1
                    DataInt(i, j, k) = Data(i, j, k)
                Next
            Next
        Next
        Refreshdisplay()
    End Sub
    Private Sub OK_Button_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles OK_Button.Click
        Select Case CheckedListBox1.SelectedIndex
            Case 0
                '              Me.DialogResult = System.Windows.Forms.DialogResult.OK
                '              Dim a As New Threading.Thread(AddressOf 
                Me.ComputeMeanAndCorrect()
                '              a.Priority = Threading.ThreadPriority.Highest
                '              a.Start()
            Case 1
                '                Me.DialogResult = System.Windows.Forms.DialogResult.OK
                '                Dim a As New Threading.Thread(AddressOf 
                Me.ComputeMeanAndCorrect2()
                '               a.Priority = Threading.ThreadPriority.Highest
                '              a.Start()
        End Select
        
    End Sub

    Sub Refreshdisplay(Optional ByVal slice As Integer = 0)
        If slice = 0 Then slice = CInt(Math.Round(DataInt.z / 2))
        CustomPictureBox1.Image = DataInt.SliceToBitmap(DataInt.TransverseSlice(slice))
        CustomPictureBox1.Refresh()
    End Sub
    Private Sub Cancel_Button_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Cancel_Button.Click
        Me.DialogResult = System.Windows.Forms.DialogResult.Cancel
        If MsgBox("Replace Anatomy data permanently?", MsgBoxStyle.YesNo, "Close") = MsgBoxResult.Yes Then
            FormMain.AnatomyVolume.Matrix = Me.DataInt
        End If
        FormMain.AnatomyDisplay.OpenFile(FormMain.AnatomyVolume, "Anatomy", True, True)
        Me.Close()
    End Sub

    Private WindowSize As Integer
    Private DataInt As MatrixDataStructures.Matrix3DInt16
    Private DataSng As MatrixDataStructures.Matrix3DSingle
    
    Sub ComputeMeanAndCorrect()
        'Calculate GlobalMean & MinThrehsold
        Dim CutOffPercentage As Integer
        If Int32.TryParse(TextBox2.Text, WindowSize) = True Then WindowSize = CInt(TextBox2.Text) Else Exit Sub
        If Int32.TryParse(TextBox1.Text, CutOffPercentage) = False Then Exit Sub
        Dim MinThreshold As Double = ((DataInt.Maximum - DataInt.Minimum) * CutOffPercentage / 100) + DataInt.Minimum
        Dim GlobalMean As Long = 0
        Dim SkippedCount As Integer = 0
        For i As Integer = 0 To DataInt.x - 1
            For j As Integer = 0 To DataInt.y - 1
                For k As Integer = 0 To DataInt.z - 1
                    If DataInt(i, j, k) > MinThreshold Then GlobalMean += DataInt(i, j, k) Else SkippedCount += 1
                Next
            Next
        Next
        GlobalMean = CLng(Math.Round(GlobalMean / ((DataInt.x * DataInt.y * DataInt.z) - SkippedCount)))
        'Now Scale Data
        Dim HalfKernelSize As Integer = CInt((WindowSize - 1) / 2)
        Dim LocalMean, Count As Integer
        For i As Integer = 0 To DataInt.x - 1
            Me.SetText("Inhomogenity Correction " & " - " & CInt(i * 100 / DataInt.x) & "% Complete.")
            For j As Integer = 0 To DataInt.y - 1
                For k As Integer = 0 To DataInt.z - 1
                    If DataInt(i, j, k) < MinThreshold Then Exit For
                    'now for each point in the datamatrix
                    LocalMean = 0
                    'start fresh from the left-topmost of the kernel
                    Count = 0
                    For ki As Integer = -1 * HalfKernelSize To 1 * HalfKernelSize
                        For kj As Integer = -1 * HalfKernelSize To 1 * HalfKernelSize
                            For kk As Integer = -1 * HalfKernelSize To 1 * HalfKernelSize
                                'skip points that are outside the image when the kernel convolves 
                                'with the data matrix
                                If i + ki < 0 Or i + ki > DataInt.x - 1 Or _
                                j + kj < 0 Or j + kj > DataInt.y - 1 Or _
                                k + kk < 0 Or k + kk > DataInt.z - 1 _
                                Then
                                    Exit For
                                End If
                                If DataInt(i + ki, j + kj, k + kk) < MinThreshold Then
                                    Exit For
                                Else
                                    LocalMean += CInt(DataInt(i + ki, j + kj, k + kk))
                                    Count += 1
                                End If
                            Next
                        Next
                    Next
                    DataInt(i, j, k) *= CShort(GlobalMean / (LocalMean / Count))
                Next
            Next
            Refreshdisplay()
        Next
        Me.SetText("Inhomogenity Correction ")
        Me.Refreshdisplay()

    End Sub
    Sub ComputeMeanAndCorrect2()
        'Calculate GlobalMean & MinThrehsold
        Dim CutOffPercentage As Integer
        If Int32.TryParse(TextBox2.Text, WindowSize) = True Then WindowSize = CInt(TextBox2.Text) Else Exit Sub
        If Int32.TryParse(TextBox1.Text, CutOffPercentage) = False Then Exit Sub
        Dim MinThreshold As Double = ((DataInt.Maximum - DataInt.Minimum) * CutOffPercentage / 100) + DataInt.Minimum
        Dim GlobalMean As Long = 0
        Dim SkippedCount As Integer = 0
        '  me.settext("Inhomogenity Correction :: Creating Mask... with threshold " & MinThreshold

        Dim scaledtemp As New Matrix3DInt16(DataInt.x, DataInt.y, DataInt.z)
        Dim max As Integer = DataInt.Maximum
        Dim min As Integer = DataInt.Minimum
        For i As Integer = 0 To DataInt.x - 1
            For j As Integer = 0 To DataInt.y - 1
                For k As Integer = 0 To DataInt.z - 1
                    scaledtemp(i, j, k) = CShort((DataInt(i, j, k) - min) * 255 / (max - min))
                Next
            Next
        Next

        For i As Integer = 0 To DataInt.x - 1
            For j As Integer = 0 To DataInt.y - 1
                For k As Integer = 0 To DataInt.z - 1
                    If DataInt(i, j, k) > MinThreshold Then GlobalMean += scaledtemp(i, j, k) Else SkippedCount += 1
                Next
            Next
        Next

        GlobalMean = CLng(Math.Round(GlobalMean / ((DataInt.x * DataInt.y * DataInt.z) - SkippedCount)))

        Dim mask As New MatrixDataStructures.Matrix3DInt16(DataInt.x, DataInt.y, DataInt.z)
        For i As Integer = 0 To DataInt.x - 1
            For j As Integer = 0 To DataInt.y - 1
                For k As Integer = 0 To DataInt.z - 1
                    If DataInt(i, j, k) > MinThreshold Then mask(i, j, k) = scaledtemp(i, j, k) Else mask(i, j, k) = CShort(GlobalMean)
                Next
            Next
        Next
        CustomPictureBox1.Image = mask.SliceToBitmap(mask.TransverseSlice(CInt(mask.z / 2)))
        '     CustomPictureBox1.Refresh()

        Me.SetText("Inhomogenity Correction :: Smoothing Mask...")
        If MsgBox("Use large kernel size (Yes)? Or large iterations (No)?", MsgBoxStyle.YesNo, "Options") = MsgBoxResult.Yes Then
            Dim a As New MatrixDataStructures.Gaussian3D(WindowSize, 0.5)
            mask = ApplyGaussian(mask, CInt(WindowSize / 2), 1.5, True)
        Else
            Dim resizer As New MatrixResizer.MatrixResize(1)
            mask = resizer.Resize(mask, CInt(mask.x * 3 / 8), CInt(mask.y * 3 / 8), CInt(mask.z * 3 / 8))
            For i As Integer = 0 To WindowSize - 1
                Me.SetText("Inhomogenity Correction :: Smoothing Mask... " & i & "/" & WindowSize)
                mask = ApplyGaussian(mask, 2, 1.5, False)
                CustomPictureBox1.Image = mask.SliceToBitmap(mask.TransverseSlice(CInt(mask.z / 2)))
                '     CustomPictureBox1.Refresh()
            Next
            mask = resizer.Resize(mask, DataInt.x, DataInt.y, DataInt.z)
        End If

        ' FormMain.AnatomyDisplay.OpenFile(mask, "Inhomogenity Mask")

        Me.SetText("Inhomogenity Correction :: Applying Mask...")

        scaledtemp.Dispose()

        Dim datasng As New Matrix3DSingle(DataInt.x, DataInt.y, DataInt.z)
        For i As Integer = 0 To DataInt.x - 1
            For j As Integer = 0 To DataInt.y - 1
                For k As Integer = 0 To DataInt.z - 1
                        datasng(i, j, k) = CSng((GlobalMean * DataInt(i, j, k) / (mask(i, j, k))))
                Next
            Next
        Next
        Dim maxs As Single = datasng.Maximum
        MsgBox(maxs & "," & max)

        GlobalMean = 0


        For i As Integer = 10 To DataInt.x - 11
            SetText("Inhomogenity Correction :: Applying Mask " & i & "/" & DataInt.x)
            For j As Integer = 10 To DataInt.y - 11
                For k As Integer = 10 To DataInt.z - 11

                    DataInt(i, j, k) = Convert.ToInt16(datasng(i, j, k) * 10000 / maxs)
                    GlobalMean += DataInt(i, j, k)

                Next
            Next
        Next

        GlobalMean = CLng(GlobalMean / ((DataInt.x - 20) * (DataInt.y - 20) * (DataInt.z - 20)))
        For i As Integer = 0 To 10
            For j As Integer = 0 To 10
                For k As Integer = 0 To 10
                    DataInt(i, j, k) = CShort(GlobalMean)
                    DataInt(DataInt.x - 1 - i, DataInt.y - 1 - j, DataInt.z - 1 - k) = CShort(GlobalMean)
                Next
            Next
        Next

        
        SetText("Inhomogenity Correction :: Done.")

        Me.Refreshdisplay()

    End Sub
    Private Function CreateKernel(ByVal halfkernelsize As Integer, ByVal variance As Double) As Matrix3DSingle
        Dim Kernel As New Matrix3DSingle(halfkernelsize * 2 + 1, halfkernelsize * 2 + 1, halfkernelsize * 2 + 1)
        Dim val As Double = 0.0

        For i As Integer = -1 * halfkernelsize To 1 * halfkernelsize
            For j As Integer = -1 * halfkernelsize To 1 * halfkernelsize
                For k As Integer = -1 * halfkernelsize To 1 * halfkernelsize
                    val = (i * i + j * j + k * k) / (-2 * variance)
                    Kernel(i + halfkernelsize, j + halfkernelsize, k + halfkernelsize) = CSng(Math.Exp(val) / ((2 * Math.PI * variance) ^ 0.5))
                Next
            Next
        Next
        Return Kernel
    End Function
    Dim Kernel As Matrix3DSingle
    Public Function ApplyGaussian(ByVal DataMatrix As Matrix3DInt16, ByVal halfkernelsize As Integer, ByVal variance As Double, Optional ByVal updatewindow As Boolean = False) As Matrix3DInt16
       
        If Kernel Is Nothing Then Kernel = CreateKernel(halfkernelsize, variance)
        Dim temp As New Matrix3DInt16(DataMatrix.x, DataMatrix.y, DataMatrix.z)
        For i As Integer = 0 To DataMatrix.x - 1
            For j As Integer = 0 To DataMatrix.y - 1
                For k As Integer = 0 To DataMatrix.z - 1
                    If DataMatrix.data(i, j, k) > 0 Then
                        'now for each point in the datamatrix
                        'start fresh from the left-topmost of the kernel
                        For ki As Integer = -1 * halfkernelsize To 1 * halfkernelsize
                            For kj As Integer = -1 * halfkernelsize To 1 * halfkernelsize
                                For kk As Integer = -1 * halfkernelsize To 1 * halfkernelsize
                                    'skip points that are outside the image when the kernel convolves 
                                    'with the data matrix
                                    If i + ki < 0 Or i + ki > DataMatrix.x - 1 Or _
                                    j + kj < 0 Or j + kj > DataMatrix.y - 1 Or _
                                    k + kk < 0 Or k + kk > DataMatrix.z - 1 Then
                                        'tempvalue = 0
                                        Exit For
                                    End If
                                    temp(i, j, k) = CShort(temp(i, j, k) + DataMatrix.data(i + ki, j + kj, k + kk) * Kernel.data(ki + halfkernelsize, kj + halfkernelsize, kk + halfkernelsize))
                                Next
                            Next
                        Next
                    End If
                Next
            Next
        Next
        Dim scaledtemp As New Matrix3DInt16(DataMatrix.x, DataMatrix.y, DataMatrix.z)
        Dim max As Integer = temp.Maximum
        Dim min As Integer = temp.Minimum
        For i As Integer = 0 To DataMatrix.x - 1
            For j As Integer = 0 To DataMatrix.y - 1
                For k As Integer = 0 To DataMatrix.z - 1
                    scaledtemp(i, j, k) = CShort((temp(i, j, k) - min) * 254 / (max - min) + 1)
                Next
            Next
        Next

        Return scaledtemp
    End Function

    Delegate Sub TextCallBack(ByVal [text] As String)

    Private Sub SetText(ByVal [text] As String)
        If Me.InvokeRequired = True Then
            Dim d As New TextCallBack(AddressOf SetText)
            Me.Invoke(d, New Object() {[text]})
        Else
            Me.Text = [text]
        End If
    End Sub

End Class
